Codecov Report

Attention: Patch coverage is 99.32432% with 1 line in your changes missing coverage. Please review.

Project coverage is 98.87%. Comparing base (c12d7bc) to head (990f9c5).

@@            Coverage Diff             @@
##             main     #174      +/-   ##
==========================================
+ Coverage   98.85%   98.87%   +0.01%     
==========================================
  Files          25       25              
  Lines        4180     4160      -20     
==========================================
- Hits         4132     4113      -19     
+ Misses         48       47       -1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

Heads up, DI v0.6.46 supports nested tuples and named tuples of arrays as caches

Thanks for the info, right now I do:

mycontext = (
    Constant(a),
    Cache(b),
    Cache(c),
    Cache(d),
)
grad_prep = prepare_gardient!(myfunc, grad_backend, x, mycontext...; strict=Val(true))

and a similar splatting for the context in gradient!. Do you this splatting is problematic for the performance ? I did some quick benchamarks and I measured no impact at all.

I don't think it is problematic, but for some users it is a bit painful: when you have tens of caches you may want to structure the more neatly, so I give that option too.

Alright thanks! I thinks it's neat enough in the code like that.

@baggepinnen It's ready to be merged, if you have some comments. You can have a look a NonLinMPC and MovingHorizonEstimator documentation (additionnal details in their respective Extended Help)

edit: I'm running the benchmarks of our papers in this version, and there is no significative difference on the performances, good news ❤️

edit: I'm running the benchmarks of our papers in this version, and there is no significative difference on the performances, good news ❤️

If you want to play around some more, a parameter that might impact performance is the coloring algorithm used by SparseMatrixColorings, more specifically the vertex ordering. Feel free to switch the default NaturalOrder() and see if it makes a difference with SMC.ncolors(jacobian_prep). Roughly speaking, the sparse autodiff complexity scales with the number of colors $c$ (or with $\lceil c / 12 \rceil$ for ForwardDiff).

If you want to play around some more, a parameter that might impact performance is the coloring algorithm used by SparseMatrixColorings, more specifically the vertex ordering. Feel free to switch the default NaturalOrder() and see if it makes a difference with SMC.ncolors(jacobian_prep). Roughly speaking, the sparse autodiff complexity scales with the number of colors c (or with ⌈ c / 12 ⌉ for ForwardDiff).

I did not played a lot with sparse jacobians for now. The pendulum exemples in our paper was small in term of number of decision variables, thus it is expected that no matter the chosen coloring/compression/etc algorithm, it will always be a bit slower then a dense formulation (for the same AD backend). It's like classical matrix operations e.g. product: if the matrices are relatively small, dense operations are generally faster even if they include many zeros.

And that was indeed the case from my quick tests. That is, on the inverted pendulum model with a control horizon of $H_c=2$ and a prediction horizon of $H_p=20$, the fastest to the slowest:

1. Single shooting transcription (dense formulation of the problem) with dense ForwardDiff, about 12 ms per time step
2. Multiple shooting transcription (sparse formulation of the problem) with dense ForwardDiff, about 20 ms per time step
3. Multiple shooting transcription (sparse formulation of the problem) with sparse ForwardDiff, about 29 ms per time step

I may improve a little bit 3. by tweaking the colouring algorithm here, but IMO that's not the good case study to adequately benchmark sparse computations.

@gdalle I just read that for Constant() context some backends requires Number or AbstractArray.

Right now I pass a mpc::NonLinMPC object in the Constant arg. I understand that it won't work with some backends ? Do you know which one ?

I think it will only be symbolic backends. I should update the docstring. Not a hard limitation either, I just didn't get around to it.

Okay, I just tested it and yes AutoSymbolics does not work with non-AbstractArray in Constant. But it does not support Cache neither, and I use them extensively also (with normal Vectors).

I can use a closure and capture mpc to avoid passing it as a Constant. There is a small impact on the performance tho. But I don't have a workaround my Cached vectors (except by manually reallocating at each call).

I will leave it as is for now.

But it does not support Cache neither, and I use them extensively also (with normal Vectors).

Did you find a bug there? I think the docs might be out of date, Caches with numbers or arrays should be supported by Symbolics too:

julia> using DifferentiationInterface

julia> import Symbolics

julia> f(x, c) = sum(copyto!(c, x))
f (generic function with 1 method)

julia> gradient(f, AutoSymbolics(), [1.0], Cache([0.0]))
1-element Vector{Int64}:
 1

Okay, I just tested it and yes AutoSymbolics does not work with non-AbstractArray in Constant.

The reason why symbolic backends (and ReverseDiff) are more difficult for handling constants is that they generate an executable (or tape) for the function during preparation. If we close over the constants during preparation, they will be hardcoded into the executable. Therefore, we need to turn the constants into symbolic variables as well, so that the user is free to give different values at runtime.

I guess theoretically we could have a separation between "constant" Constants (those which are the same at preparation and execution) and "variable" Constants (those which differ between preparation and execution). But you can probably see how that would be confusing for the casual user, unless we come up with better names ^^

But it does not support Cache neither, and I use them extensively also (with normal Vectors).

Did you find a bug there? I think the docs might be out of date, Caches with numbers or arrays should be supported by Symbolics too:

julia> using DifferentiationInterface

julia> import Symbolics

julia> f(x, c) = sum(copyto!(c, x))
f (generic function with 1 method)

julia> gradient(f, AutoSymbolics(), [1.0], Cache([0.0]))
1-element Vector{Int64}:
 1

Edit: I remember now, I marked Caches as non-supported by Symbolics because tests fail for the Jacobian only:

https://github.com/JuliaDiff/DifferentiationInterface.jl/blob/cfab84dd775f6b7e9b39d0d0f84f0d260fc4fe09/DifferentiationInterface/test/Back/SymbolicBackends/symbolics.jl#L22-L27

But it does not support Cache neither, and I use them extensively also (with normal Vectors).

Did you find a bug there? I think the docs might be out of date, Caches with numbers or arrays should be supported by Symbolics too:

julia> using DifferentiationInterface

julia> import Symbolics

julia> f(x, c) = sum(copyto!(c, x))
f (generic function with 1 method)

julia> gradient(f, AutoSymbolics(), [1.0], Cache([0.0]))
1-element Vector{Int64}:
 1

Maybe, AutoSymbolics simply throws an error is I leave the Constant(mpc). If I remove the Constant(mpc) but leave the multiples Cache(vector_a), Cache(vector_b), ... it return a derivative but it seems imprecise since Ipopt suddenly returns error code like "NUMERICAL_ERROR" and "ALMOST_LOCALLY_SOLVED" stuff like that (or maybe it is precise at some points of decision variable and not other, hard to say) . There is no error code at all with AutoForwardDiff and AutoFiniteDiff.

Hard to troubleshoot this kind of problem. It feels like I just entered a rabit hole haha.

Okay, I just tested it and yes AutoSymbolics does not work with non-AbstractArray in Constant.

The reason why symbolic backends (and ReverseDiff) are more difficult for handling constants is that they generate an executable (or tape) for the function during preparation. If we close over the constants during preparation, they will be hardcoded into the executable. Therefore, we need to turn the constants into symbolic variables as well, so that the user is free to give different values at runtime.

I guess theoretically we could have a separation between "constant" Constants (those which are the same at preparation and execution) and "variable" Constants (those which differ between preparation and execution). But you can probably see how that would be confusing for the casual user, unless we come up with better names ^^

Yeah, this "nomenclature issue" crossed my mind during the implementation. And I literally use both "defintion" in the package right now: Constant(mpc) is literally a constant that won't change at runtime.

In the linearize implementation of #179, the Constants are vectors that are not mutated during the differentiation but their value is different at each jacobian! call.

Worst idea ever: a new context called ConstantConstant() 😆

Okay, I'm a bit ambivalent on if should keep the Constant(mpc) or not.

My feeling right now: since it does not need any modification at runtime (same object at each jacobian! call), it's a "constant" Constant and I can use the closure trick. Theoretically, if AutoSymbolics improve in the future, Cache with normal vectors will be well supported, and it will work well for gradient, jacobian etc.

Conclusion: I will remove Constant(mpc) (and same thing for MHE).

My feeling right now: since it does not need any modification at runtime (same object at each jacobian! call), it's a "constant" Constant and I can use the closure trick.

I tend to agree with that. Essentially, if your constant never changes, it's part of your function, and then why should DI give a shit ^^

Theoretically, if AutoSymbolics improve in the future, Cache with normal vectors will be well supported, and it will work well for gradient, jacobian etc.

Looking at it again, I think it is a bug in Symbolics. Will try to dig up an MWE

@baggepinnen I will merge the PR, you can comment here if you spot anything not clear or sus